rm(list = ls())
# Function to load packages
loadPkg=function(toLoad){
       for(lib in toLoad){
              if(! lib %in% installed.packages()[,1])
              {install.packages(lib, repos='http://cran.rstudio.com/')}
              suppressMessages( library(lib, character.only=TRUE))}}

# Load libraries
packs=c("maps", "rworldmap", "maptools", "sp", "spdep", "gstat", "cshapes","readstata13", "data.table",
        "spatstat", "DataCombine", "countrycode",  "cshapes", "arm", "rgeos", "raster", "gpclib", "tidyr", "pwt9",
        "texreg", "reshape2","ggplot2", 'foreign', 'car', 'lme4', 'dplyr', "mapdata", "rgdal","gridExtra",
        "pscl", "stargazer","stringi","multiwayvcov", "lmtest","scales","stringi","haven",
        "coefplot", "ggthemes", "gridExtra", "grid", "cowplot", "GGally")
loadPkg(packs)

##########################################
##########################################
setwd("YOURPATH")
dir.create("figures")
##########################################
##########################################

### I make a set of R functions for plotting. To replicate the all figures,
#you need to load them to make the plots with a simulation approach 
source("setup_code.R")

#set ggplot theme
theme_set(theme_gray())
#load the data
load("paperdata.RData")


data <- paperdata %>% 
       dplyr::filter(weaker_chal==1)

###########################################################################
############### This is the main results in the analysis (Models 1 - 6)
###########################################################################
#Model 1: base
b1 <- bayesglm(initiation ~ tgt_landneighbors2 +  backdown_terr1 +  backdown_nonterr1 + 
                      joint_demo +  defensetgt + offensechal +  land_neighbor + 
                      vicratio +  tmid +  tmid_sq + tmid_cub, 
          data = data, family = binomial(link = "logit"),
          prior.scale=2.5, prior.df=1)

# model 2: 
b2 <- bayesglm(initiation ~ tgt_landneighbors2 +  backdown_terr1 +  backdown_nonterr1 + 
                      joint_demo +  pot_milicapchal +  pot_milicaptgt +  land_neighbor + 
                     tmid +  tmid_sq + tmid_cub, 
               data = data, family = binomial(link = "logit"),
               prior.scale=2.5, prior.df=1)

#model 3
b3 <- bayesglm(initiation ~ tgt_landneighbors2 +  backdown_terr1 +  backdown_nonterr1 + 
                      joint_demo +  intra_wartgt +  inter_war_totaltgt + defensetgt +
                      offensechal +  land_neighbor + 
                      vicratio +  tmid +  tmid_sq + tmid_cub, 
               data = data, family = binomial(link = "logit"),
               prior.scale=2.5, prior.df=1)
# Model 4: 
b4 <- bayesglm(initiation ~ tgt_landneighbors2 +  backdown_terr1 +  backdown_nonterr1 + 
                      joint_demo +  solschangetgtdum + solschangechaldum +
                      offensechal +  land_neighbor + defensetgt  +
                      vicratio +  tmid +  tmid_sq + tmid_cub, 
               data = data, family = binomial(link = "logit"),
               prior.scale=2.5, prior.df=1)

# Model 5 
b5 <- bayesglm(terrmidl ~ tgt_landneighbors2 +  backdown_terr1 +  backdown_nonterr1 + 
                      joint_demo +  offensechal +  land_neighbor + defensetgt  +
                      vicratio +  tterr + tterr2 + tterr3, 
               data = data, family = binomial(link = "logit"),
               prior.scale=2.5, prior.df=1)

#Model 6:
b6 <- bayesglm(initiation ~ tgt_landneighbors2 +  backdown_terr1 +  backdown_nonterr1 + 
                      joint_demo +  offensechal +  land_neighbor + defensetgt  +
                      weaker_chal +  major_major +  major_minor +  minor_major + 
                      vicratio +   tmid +  tmid_sq + tmid_cub, 
               data = paperdata, family = binomial(link = "logit"),
               prior.scale=2.5, prior.df=1)

################### Making coefficient plot for each model################### 
################### Figure 2 (a)-(f)
################### 
###make coefficient plot for all models

################### Figure 2(a)###################
coef_f1 <- coef_clusterplot(ModelResults = b1,
                            varnames = c("tgt_landneighbors2","backdown_terr1",
                                         "backdown_nonterr1", "joint_demo",
                                         "defensetgt", "offensechal",
                                         "land_neighbor","vicratio"),
                            data = data, 
                            clusterid = "ddyad") + 
       scale_x_discrete(labels=c("Past yielding (nonterritory)", "Past yielding",
                                 "Defensive alliance (target)", "Joint democracy",
                                 "Land neighbor", "Offensive alliance (challenger)",
                                 "Num. of land neighbors",
                                 "Exp. victory prob. ratio")) + 
       theme(axis.text = element_text(size=14),
             plot.title = element_text(hjust = 0.5,size = 30, face = "bold"),
             axis.title=element_text(size=14)) + 
       ggtitle("") 

ggsave("./figures/cof1.pdf", plot = coef_f1,  width = 10, height = 8)

################### Figure 2(b)###################
coef_f2 <- coef_clusterplot(ModelResults = b2,
                            varnames = c("tgt_landneighbors2","backdown_terr1",
                                         "backdown_nonterr1", "joint_demo",
                                         "pot_milicapchal", "pot_milicaptgt",
                                         "land_neighbor"),
                            data = data, 
                            clusterid = "ddyad") + 
       scale_x_discrete(labels=c("Past yielding (nonterritory)", "Past yielding",
                                 "Joint democracy",
                                 "Land neighbor", "Alliance strength (challenger)",
                                 "Alliance strength (target)", 
                                 "Num. of land neighbors")) + 
       theme(axis.text = element_text(size=14),
             plot.title = element_text(hjust = 0.5,size = 30, face = "bold"),
             axis.title=element_text(size=14)) + 
       ggtitle("") 
ggsave("./figures/cof2.pdf", plot = coef_f2,  width = 10, height = 8)

################### Figure 2(c)###################
coef_f3 <- coef_clusterplot(ModelResults = b3,
                            varnames = c("tgt_landneighbors2","backdown_terr1",
                                         "backdown_nonterr1", "joint_demo",
                                         "defensetgt", "offensechal",
                                         "intra_wartgt", "inter_war_totaltgt",
                                         "land_neighbor","vicratio"),
                            data = data, 
                            clusterid = "ddyad") + 
       scale_x_discrete(labels=c("Past yielding (nonterritory)", "Past yielding",
                                 "Defensive alliance (target)","Interstate war (target)",
                                 "Civil war (target)", "Joint democracy",
                                 "Land neighbor", "Offensive alliance (challenger)",
                                 "Num. of land neighbors",
                                 "Exp. victory prob. ratio")) + 
       theme(axis.text = element_text(size=14),
             plot.title = element_text(hjust = 0.5,size = 30, face = "bold"),
             axis.title=element_text(size=14)) + 
       ggtitle("") 

ggsave("./figures/cof3.pdf", plot = coef_f3,  width = 10, height = 8)

################### Figure 2(d) ###################
coef_f4 <- coef_clusterplot(ModelResults = b4,
                            varnames = c("tgt_landneighbors2","backdown_terr1",
                                         "backdown_nonterr1", "joint_demo",
                                         "defensetgt", "offensechal",
                                         "solschangetgtdum", "solschangechaldum",
                                         "land_neighbor","vicratio"),
                            data = data, 
                            clusterid = "ddyad") + 
       scale_x_discrete(labels=c("Past yielding (nonterritory)", "Past yielding",
                                 "Defensive alliance (target)",
                                 "Joint democracy",
                                 "Land neighbor", "Offensive alliance (challenger)",
                                 "SOLS change (challenger)",
                                 "SOLS change (target)",
                                 "Num. of land neighbors",
                                 "Exp. victory prob. ratio")) + 
       theme(axis.text = element_text(size=14),
             plot.title = element_text(hjust = 0.5,size = 30, face = "bold"),
             axis.title=element_text(size=14)) + 
       ggtitle("") 

ggsave("./figures/cof4.pdf", plot = coef_f4,  width = 10, height = 8)

################### Figure 2(e)###################
coef_f5 <- coef_clusterplot(ModelResults = b5,
                            varnames = c("tgt_landneighbors2","backdown_terr1",
                                         "backdown_nonterr1", "joint_demo",
                                         "defensetgt", "offensechal",
                                         "land_neighbor","vicratio"),
                            data = data, 
                            clusterid = "ddyad") + 
       scale_x_discrete(labels=c("Past yielding (nonterritory)", "Past yielding",
                                 "Defensive alliance (target)", "Joint democracy",
                                 "Land neighbor", "Offensive alliance (challenger)",
                                 "Num. of land neighbors",
                                 "Exp. victory prob. ratio")) + 
       theme(axis.text = element_text(size=14),
             plot.title = element_text(hjust = 0.5,size = 30, face = "bold"),
             axis.title=element_text(size=14)) + 
       ggtitle("") 

ggsave("./figures/cof5.pdf", plot = coef_f5,  width = 10, height = 8)

################### Figure 2(f)###################
coef_f6 <- coef_clusterplot(ModelResults = b6,
                            varnames = c("tgt_landneighbors2","backdown_terr1",
                                         "backdown_nonterr1", "joint_demo",
                                         "defensetgt", "offensechal",
                                         "land_neighbor","vicratio",
                                         "weaker_chal", "major_major", "major_minor",
                                         "minor_major"),
                            data = paperdata, 
                            clusterid = "ddyad") + 
       scale_x_discrete(labels=c("Past yielding (nonterritory)", "Past yielding",
                                 "Defensive alliance (target)", "Joint democracy",
                                 "Land neighbor","Major--Major powers","Major--Minor powers",
                                 "Minor--Major powers", 
                                 "Offensive alliance (challenger)",
                                 "Num. of land neighbors",
                                 "Exp. victory prob. ratio",
                                 "Weaker challenger dyad")) + 
       theme(axis.text = element_text(size=14),
             plot.title = element_text(hjust = 0.5,size = 30, face = "bold"),
             axis.title=element_text(size=14)) + 
       ggtitle("") 

ggsave("./figures/cof6.pdf", plot = coef_f6,  width = 10, height = 8)



################### Make an ROC plot ################### 
################### Figure 6 ########################
####################################################


b0 <- bayesglm(initiation ~ joint_demo +  defensetgt + offensechal +  
                      vicratio +  tmid +  tmid_sq + tmid_cub, 
               data = data, family = binomial(link = "logit"),
               prior.scale=2.5, prior.df=1)
roc <- plot_roc(ModelResults = list(b0, b1), 
                linetypes = c("solid", "dotted"), interval=0.2)
ggsave("./figures/roc.png", plot = roc)


